using System;
using System.IO;

namespace NumericalSolution {
  class ErrorComparator {

    static string outfile = "c:/nsworkspace/diffusion1d_imp_error.csv";

    public static void Main() {

      int initialNumberOfNodes = 11;
      int endNumberOfNodes = 81;
      int stepNumberOfNodes = 5;

      double initialDeltaTime = 0.001;
      double endDeltaTime = 0.005;
      double stepDeltaTime = 0.001;

      StreamWriter writer = File.CreateText(outfile);

      writer.Write("#Nodes, ");
      for (double dt = initialDeltaTime; 
           dt <= endDeltaTime; 
           dt += stepDeltaTime) {
        writer.Write("dt(" + dt + "), ");
      }
      writer.WriteLine();

      for (int n = initialNumberOfNodes; 
           n <= endNumberOfNodes; 
           n += stepNumberOfNodes) {

        writer.Write(n + ", ");
        
        for (double dt = initialDeltaTime; 
             dt <= endDeltaTime; 
             dt += stepDeltaTime) {

          Diffusion1D exact = new Diffusion1DExact();
          Diffusion1D numerical = new Diffusion1DImplicit();

          exact.NumberOfNodes = n;
          numerical.NumberOfNodes = n;
          exact.DeltaTime = dt;
          numerical.DeltaTime = dt;
          exact.EndTime = 1.0;
          numerical.EndTime = 1.0;

          double error = CalculateError(exact, numerical);
          writer.Write(error + ", ");
        }
        writer.WriteLine();

      }
      writer.Close();
    }

    static private double CalculateError(Diffusion1D exact, 
                                         Diffusion1D numerical) {

      double[] u_exact;
      double[] u_numerical;
      double error = 0.0;

      int count = 0;
      for (exact.CurrentTime = 0.0;
           exact.CurrentTime <= exact.EndTime;
           exact.CurrentTime += exact.DeltaTime) {

        if (exact.CurrentTime == 0.0) {
          u_exact = exact.Initialize();
          u_numerical = numerical.Initialize();
          continue;
        }

        u_exact = exact.Next();
        u_numerical = numerical.Next();

        for (int i = 1; i < u_exact.Length - 1; i++) {
          error += Math.Abs(u_exact[i] - u_numerical[i]) / u_exact[i];
          count++;
        }
      }
      return error / (double)count;
    }
  }
}
